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ABSTRACT 

Differential GPS (DGPS) positioning is used to accurately locate a GPS receiver based 
upon the well-known position of a reference site. In utilizing this technique, several 
errors sources contribute to position inaccuracy. This thesis investigates the error in 
DGPS operation and attempts to develop a statistical model for the behavior of this 
error. The model for DGPS error is developed using GPS data collected by Draper 
Laboratory. The Marquardt method for non-linear curve-fitting is used to find the 
parameters of a first order Markov process that models the average errors from the col- 
lected data. The results show that a first order Markov process can be used to model 
the DGPS error as a function of baseline distance and time delay. The model’s time* 
correlation constant is 3847.1 seconds (1.07 hours) for the mean square error. The dis- 
tance correlation constant is 122.8 kilometers. The total process variance for the 
DGPS model is 3.73 meters^. 
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Chapter 1 
Introduction 

1.1 Motivation 

The Global Positioning System (GPS) is used for navigation and associated tasks 
in a wide variety of applications. One of the primary modes of GPS operation is differ- 
ential GPS. Differential GPS allows more accurate positioning than stand-alone GPS. 
Many new applications are possible because of the increased accuracy obtained 
through differential GPS. For instance, a cruise missile operating with differential GPS 
can be guided very accurately to its target destination. 

The benefits of differential GPS are apparent, but it must be well understood in 
order to maximize this benefit. There is a need, therefore, to determine the positioning 
errors that are associated with differential GPS operation for a variety of cases, then 
the user can make knowledgeable decisions on the implementation of the differential 
technique. 

1.2 GPS Description 

GPS signals are transmitted from the satellites at two different frequencies, LI at 
1575.42 MHz and L2 at 1227.60 MHz. The LI frequency is modulated by two pseu- 
dorandom noise codes, C/A-code at 1.023 MHz and P-code at 10.23 MHz. The L2 car- 
rier is only modulated by the P-code. Since the wavelength of the P-code is 30 meters 
compared to 300 meters for the C/A-code, the P-code provides much better accuracy 
for positioning. The P-code became known as the Y-code when it was encrypted so 
only authorized users could access it. A P(Y)-code capable receiver must be “keyed” 
with the proper information to decode the incoming signals and obtain P-code accu- 
racy. [7] 



For C/A“Code users, the performance is degraded due to a purposeful degradation 
in the quality of broadcast orbits and satellite clock dithering. This degradation is 
called selective availability (SA). [24] An advantage of P-code operation is that per- 
formance is not affected by selective availability. 

Positioning is accomplished by determining the time it takes for the GPS signals 
from each satellite to reach the receiver. This time provides an estimate of the distance 
from satellite to receiver, and “triangulating” using the signals from four or more satel- 
lites provides the receiver’s position in three dimensions. The transmit and receive 
times are corrupted by clock errors, so a true range measurement is not possible. 
Therefore, the distance measurement from the receiver to each satellite is called the 
pseudorange. 

The “triangulation” is not perfect, because the satellites are not located optimally 
with respect to the receiver. Additional error is contained in the navigation solution 
due to the position of the satellites. The instantaneous measure of this error is called 
the Geometric Dilution of Precision (GDOP). It represents the rms error realized in the 
navigation solution for unit error in the triangulation measurements. GDOP accounts 
for dilution of precision in three dimensions and time. GDOP can be broken down into 
the following components: vertical (VDOP), horizontal (HDOP), time (TDOP), posi- 
tion (PDOP). PDOP reflects the dilution of precision in three dimensions. [20] This 
can be further decomposed into XDOP, YDOP, and ZDOP. 

Incoming GPS signals are delayed when they pass through the earth’s atmosphere. 
The ionosphere is one of the primary causes for this delay. A P-code receiver can track 
both the LI and L2 signals. Since these signals are at different frequencies, the iono- 
spheric delays will be different for each. An estimate for the ionospheric delay can be 
obtained utilizing the difference between the LI and L2 delays. Thus a P-code receiver 
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provides a dual-frequency ionospheric correction term that is subtracted from pseudo- 
ranges. Since the LI and L2 delays are obtained from tracking loops on the signals, 
additional noise enters the navigation solution when the dual-frequency correction is 
used. [10] 

GPS satellites broadcast almanac files that contain the orbit parameters required to 
predict the satellites’ positions. These files are updated every week in order to main- 
tain the accuracy. All GPS data is referred to by week. The system clock is reset to 
zero every week on midnight Saturday, so the correct date is found by referencing the 
GPS week in the data files. 

1.3 Differential GPS 

As shown in Figure 1.1, positioning using differential GPS (DGPS) involves the 
use of two GPS receivers. One receiver is located at a fixed position that is accurately 
surveyed, while the other is allowed to move around. The goal is to accurately deter- 
mine the position of the “roaming” receiver. 
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Since the position of receiver A is surveyed, corrections can be determined for the 
pseudoranges that generate the GPS position. These corrections are then sent to 
receiver B. The roaming receiver uses the corrections from receiver A to update or 
correct its position. This procedure is highly accurate with small baseline lengths. 
Errors increase as the baseline increases or as a time delay is introduced between the 
measurement and use of the corrections. Such a time delay can be the result of the 
desire for covert operation. 

Using DGPS can help to eliminate errors due to the satellite clock drift, ephemeris 
errors, and propagation delays. Errors remaining include receiver interchannel biases, 
receiver noise and quantization errors, multipath, residual propagation delays, and 
errors due to time delays in the pseudorange corrections. [4] 

1.4 Previous Work 

Several sources were found that discussed DGPS errors. Most of the previous 
work was performed with C/A-code receivers. Error budgets are cited in [4, 15, 17] 
and depend on the type of receiver being used. For the most part, these budgets were 
the specifications for the devices or analytic estimates which do not always correlate 
with actual performance. The next section provides one of these error budgets. An 
analytic development of the individual error sources was provided in [9]. 

Two references were found that model GPS and DGPS errors as first order Markov 
processes. [7, 23] A Markov model was chosen because many physical phenomena 
exhibit behavior that fit such models reasonably well. No GPS data was used in deriv- 
ing these models. 

One reference used simulated data because receiver data for the desired scenarios 
was unavailable. [14] Some studies used C/A-code data for their analysis, including 
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data found on the Internet. [5, 13] One study did use P-code receivers, but this study 
examined the DGPS errors for a specific geometry and did not develop a model. [25] 
This thesis utilizes data segments with common satellites and processes them inde- 
pendently. This prevents jumps in the solutions due to satellite switches from affecting 
the results. A statistical model for the DGPS errors is provided. 

1.5 Error Budgets and Modeling 

Errors for a P-code receiver are typically cited at about 10 meters RMS in each 

axis. Table 1.1 shows the error budgets for static GPS operation, and Table 1.2 shows 
an error budget for Differential GPS. These budgets are specifications for performance 
for co-altitude receivers at a separation of 50 nmi (92.6 km). They were published by 
the Air Force GPS Range Applications Joint Program Office. [4, 15] 


Table 1.1: GPS Static Error Budget 


Error Source 

C/A-Code Predicted Error 
(meters) 

P-Code Predicted Error 
(meters) 

Satellite Clock Error 

3.05 

3.05 

Ephemeris Error 

2.62 

2.62 

Ionospheric Delay Error 

6.40 

0.40 

Tropospheric Delay Error 

.40 

0.40 

Receiver Noise 

2.44 

0.24 

Receiver Interchannel Bias 

0.61 

0.15 

Multipath 

3.05 

1.22 

UERE (RMS) 

8.54 

4.25 

RMS Horizontal Position 
Error (HDOP = 1 .5) 

12.81 

6.37 

RMS Vertical Position 
Error (VDOP = 2.5) 

21.34 

10.62 

Total RMS Error 

24.89 

12.38 
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Table 1.2: Differential GPS Error Budget 


Error Source 

C/A-Code Predicted Error 
(meters) 

P-Code Predicted Error 
(meters) 

Satellite Clock Error 

0 

0 

Ephemeris Error 

0 

0 

Residual lonospheric/Tro- 
pospheric Delay Error 

0.15 

0.15 

Receiver Noise 

2.44 

0.24 

Receiver Interchannel Bias 

0.61 

0.15 

Multipath 

3.05 

1.22 

UERE (RMS) 

3.95 

1.26 

RMS Horizontal Position 
Error (HDOP = 1.5) 

5.93 

1.89 

RMS Vertical Position 
Error (VDOP = 2.5) 

9.88 

3.15 

Total RMS Error 

11.53 

3.68 


The DGPS errors will change as a function of baseline distance and the time delay 
in the correction inputs. Previous work has estimated the change with baseline dis- 
tance to be a linear relationship. [13] One of the major challenges in this project will 
be estimating long term behavior from many segments of short term data. 

1.6 GPS Data 

Several data sets were utilized in this thesis. Draper Laboratory collected all the 
data. Three data collections were performed: East Coast experiments, West Coast 
experiments, and zero baseline, rooftop experiments. The first two experiments were 
performed under a previous effort. They provide GPS data at locations with varying 
separations so several baselines are available. The experiment names refer to the 
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region where the experiments were performed, and the objectives for both were the 
same. The rooftop experiment was performed as part of this thesis. The primary bene- 
fit was to be a characterization of the receiver noise. An unexpected benefit was char- 
acterization of errors due to multipath. This is possible because the baseline is zero, 
and both receivers take inputs from the same antenna. These experiments will be 
described in more detail in Chapter 2. The results from all three experiments will be 
combined to form a DGPS model from many different cases. 

1.7 Thesis Objectives 

This thesis will examine the errors in differential GPS as a function of baseline dis- 
tance between the two receivers and the time delay between when the corrections are 
computed and when they are used. The final results will include data from P-code 
receivers only. A model for the errors will be developed and analyzed to determine 
how well it fits the errors observed in actual differential GPS experiments. 

The thesis is organized as follows. Chapter 2 discusses the experiments and the 
data obtained. Chapter 3 presents methods used in the processing of the data and mod- 
eling of the errors. In Chapter 4, the results from each experiment is shown including 
the performance of the models developed. Chapter 5 is a summary of the thesis and 
conclusions drawn from the results. Suggestions for future study are also discussed. 
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Chapter 2 


Data Collection and Formats 

This chapter describes the data that was used to develop the model for DGPS 
errors. All of this data was collected by Draper Laboratory. Much of it comes from 
keyed P-code receivers that is unaffected by selective availability. 

2.1 Draper Rooftop Experiments 

2.1.1 Experiment Description 

The goal of the rooftop experiments was to obtain data with a zero baseline and to 
get an estimate of receiver noise. A GPS antenna was placed on a tripod on the rooftop 
of Draper. Two Trimble TANS receivers in a laboratory on the 4th floor received the 
inputs from the single antenna. Since the receivers were getting the exact same inputs, 
the only difference in their solutions was due to the receivers’ internal processing. The 
baseline is zero meters, and the time delay errors were examined independently of 
baseline dependent errors. Data was recorded every 10 seconds. 

The TANS receivers are keyed P-code receivers which are resistant to selective 
availability. The position errors are expected to be well within the specified value of 1 6 
meters. A bug in the translation program written for the TANS receivers caused errors 
in the recorded ephemeris data. A Novatel receiver was positioned on the rooftop to 
collect data for the same periods as the TANS receivers. The Novatel translation pro- 
gram worked, so it was used as the source of broadcast ephemeris data. 



2.1.2 Data Description 

The first set of data was taken on Monday, November 14th, 1995, 23:44:50 to 
Tuesday, November 15th, 22:33:10. This was GPS week 775. This set contains 
approximately 24 hours of data. Unfortunately, one of the computers controlling the 
TANS receivers failed to store the collected data. All data from one receiver was lost. 
It was thus not possible to estimate receiver noise from this data. However, the data 
from the remaining receiver is very useful for the time delay error analysis, because 
the data can be copied and shifted in time. In many trials comparing receiver A vs. 
receiver B (two different receivers) with receiver A vs. receiver A (the same receiver), 
the time shifted errors are statistically identical. Using the same receiver for the time 
delay study eliminates any sample mean in the data. 

The second set of data was taken on Thursday, November 17th, 1995, 17:02:38 to 
22:13:40. This was also GPS week 775. This set contains approximately 5 hours of 
data. This time both computers stored the data, so comparisons between the two 
receivers could be performed. Unfortunately, a gap of about 1,5 hours occurred in the 
middle of the data. The gap occurred when neither receiver was able to maintain signal 
lock on more than three satellites. The reason for the gap has not been determined. 
Possible explanations include objects blocking the line of sight to the satellite, an inac- 
tive GPS satellite, or a satellite that was not correctly transmitting at that time. It is 
unlikely, however, that any nearby structure could have blocked the satellites for 1.5 
hours. 

The receivers were configured to track the same four satellites. This was done so 
the position solution from both receivers would have common error sources. Using 
solutions from different satellites would mask the differences we are trying to observe. 
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The receivers were also configured to hold a set of four satellites until the GDOP 
became greater than 6. At that time, both receivers switched to a new set of four satel- 
lites. 


2.2 East Coast Experiments 
2.2.1 Experiment Description 

The East Coast experiments were conducted in December, 1993. These experi- 
ments included both C/A and P-code receivers. Table 2. 1 shows the location and type 
of the receivers used. Data collected from Ft. Stewart, Georgia was not used due to 

Table 2.1 East Coast Receiver Locations and Types 


Receiver Location 

Receiver 

Type 

Cambridge, MA 

P-code 

Portsmouth, NH 

P-code 

Langley AFB, VA 

P-code 

Willow Run, MI 

C/A-code 


problems with the data collection. The East Coast baselines are shown in Table 2.2. 

Table 2.2 East Coast Baselines 


Receiver A 
Location 

Receiver B 
Location 

Baseline Length 
(km) 

Cambridge, MA 

Portsmouth, NH 

84 

Langley AFB, VA 

Cambridge, MA 

740 

Langley AFB, VA 

Portsmouth, NH 

820 

Langley AFB, VA 

Willow Run, MI 

838 

Willow Run, MI 

Cambridge, MA 

1022 

Willow Run, MI 

Portsmouth, NH 

1052 
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The C/A-code receiver at Willow Run was a Novatel 95 IR. The selective availability 
in the data from this receiver was removed by Aerospace Corporation. The P-code 
receivers were Trimble 4000SSE unencrypted receivers. The P-code receivers were 
able to track P-code even though unencrypted, because P(Y)-code encryption was not 
yet in effect. 

2.2.2 Data Description 

Seventeen different sessions of data were collected. All the receivers were oper- 
ated at the same time. Some of the sessions were lost due to data collection problems. 
The data was searched for segments where all receivers had at least twenty-five min- 
utes of time tracking the same four satellites. Five such segments were found in ses- 
sions 8, 10, 14, and 16. These five data segments were used to form time delay 
differences as well as to estimate the differences due to the six baseline separations. 

2.3 West Coast Experiments 
2.3.1 Experiment Description 

The West Coast experiments were conducted in March, 1994. Keyed P-code Trim- 
ble 4000SSE receivers were used at all locations. Only two such receivers were avail- 
able, however, so each baseline had to be instrumented separately. Data was collected 
from four baselines as shown in Table 2.3. 


Table 2.3 West Coast Baselines 


Receiver A 
Location 

Receiver B 
Location 

Baseline Length 
(km) 

San Diego, CA 

Long Beach, CA 

149 

Stockton, CA 

Fallon, NV 

282 

Ft. Irwin, CA 

Fallon, NV 

496 

San Diego, CA 

Fallon, NV 

758 


24 



















2.3.2 Data Description 

Sixteen sessions of data were collected on the west coast, including four sessions 
per baseline. Data collection problems resulted in the loss of two sessions. For the 
remaining sessions, data segments of 25 minutes or greater with the two receivers 
tracking the same four satellites were selected. Each session had from 1 to 5 such seg- 
ments. These segments were again used to form the time delay differences and to esti- 
mate the error due to the four baseline separations. 

2.4 Data Anomalies 

Several interesting qualities or anomalies have been observed in the collected data 
sets. These are documented to stimulate future discussion and investigation. The gap 
in the 5-hour, zero baseline rooftop data has already been mentioned. From discus- 
sions with experienced GPS users, satellite outages have occurred on occasion. 

The GPS navigation solution jumps at every satellite switch. This is an expected 
attribute, since it is a function of all four satellites being tracked. One of the primary 
goals of the data reduction was to deal appropriately with these jumps. Figure 2. 1 
shows the data from a zero baseline 5-hour rooftop session. The graphs show the com- 
ponents of the difference between one receiver’s solution and presumed truth. The ver- 
tical lines show the places where the satellite switches occur. Notice the jumps in the 
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Z Error (meters) Y Error (meters) X Error (meters) 


output solution at the satellite breaks. Also note the large data gap in the middle of this 
set. 
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Figure 2.2 shows a close-up of the last satellite break in the zero baseline data. 
Note the fairly large jumps in both x and z directions. 




GPS Time (hours) 


Figure 2.2 Satellite Break, Zero Baseline Data 

Figure 2.3 shows a close-up of another session of zero baseline data. The vertical 
lines here are not the satellite switches. They mark each measurement made by the 
GPS receiver. These occur every 10 seconds. A satellite switch occurred at 1 14.32 sec- 
onds. Note the data gap before the switch. This occurred several times in the collected 
data, and is presumably due to a delay in locking onto the new satellite. 

The errors often exhibited cyclic behavior. A change in the frequency and magni- 
tude of these oscillations is observable at the satellite switches. The oscillations will be 
primarily attributed to multipath. The changes in the oscillations make sense since the 
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new satellite set provides a new multipath geometry. This will be discussed further in 
the next chapter. 



Figure 2.3 Satellite Switch and Data Gap, Zero Baseline Data 
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Figure 2.4 shows a sample of the West Coast data. Again the vertical lines mark 
each satellite switch. Jumps in one or more components are apparent at each satellite 
switch. 




a I u 1 1 i j I I I 

62 62.5 63 63.5 64 64.5 65 65.5 66 

GPS Time (hours) 

Figure 2.4 Sample West Coast Data 


Figure 2.5 shows a close-up of the West Coast data. The vertical lines mark each 
GPS measurement. Data was taken every 30 seconds. A satellite switch occurred at the 
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jump in the middle of the data set at 63.25 seconds. Note once again the jump in the 
solution that occurred at the satellite switch. 





GPS Time (hours) 

Figure 2.5 Ciose-up of West Coast Data 


Figure 2.6 shows an example with only the y direction plotted. Note that satellite 
switches occurred at 86.8 and 87.1 seconds. The changes in amplitude and frequency 
of the oscillations in the navigation solution are more observable in this case. Analysis 
will show that primary oscillation periods change from 7 1 .4 seconds to 90.9 seconds 
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to 125 seconds for the satellite switches in this example. This was determined by 
applying a fast Fourier transform to the data. 



Figure 2.6 Raw Data Oscillations, East Coast Data 


The GPS navigation solution was observed to jump at points other than satellite 
switches. In one case, a jump in z position of 7 meters, y of 3 meters, and x of 0.75 

meters was observed. Such jumps were rare, and the cause is unknown. One explana- 

« 

tion is that the satellite ephemeris was updated at those points. Figure 2.7 shows a 
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close-up of a jump in the solution when there was no satellite switch. The reader may 
have noticed this jump in Figure 2.1, towards the very end of the data set. 



117.6 117.605 117.61 117.615 117.62 117.625 117.63 117.635 117.64 117.645 117.65 



117.6 117.605 117.61 117.615 117.62 117.625 117.63 117.635 117.64 117.645 117.65 



GPS Time (hours) 

Figure 2.7 Navigation Solution Jump Without Satellite Switch 


Another anomaly of note was an occasional large and sometimes persistent differ- 
ence between the position solution from the two receivers. Given that the receivers 
were being driven by the same antenna, the difference had to have arisen from some 
internal difference between the two receivers. Sometimes these differences only lasted 
for one measurement, sometimes they lasted for 100 seconds (10 measurements). A 
GPS receiver is a complex non-linear device. Such behavior is not surprising and 
could arise for a number of reasons. Included are low S/N for a particular satellite or 
some receiver moding change. The level of process noise being introduced into the 
receiver Kalman filter was such that they were essentially delivering point solutions. 
This was done, of course, so that no ‘‘artificial” time correlation would be introduced 
into the data. The persistence of the anomalies cannot thus be blamed on the naviga- 
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tion filter. It should be remembered that the receivers were being forced to hold onto a 
set of satellites until the GDOP had risen to 7. This would tend to make the solution 
more sensitive to pseudo-range errors than would normally be the case. 

The ephemeris data collected by one of the Trimble receivers was missing some of 
the satellite information. It appeared to start taking data correctly, but then ephemeris 
records would show only information from about three of the satellites. Fortunately, 
the ephemeris from the Novatel receiver was available, so this loss did not affect the 
results. However, no reason for the data loss has been determined. All other elements 
of the data from that receiver seemed fine. 

The Novatel receiver ephemeris changed GPS week from 775 to 776. This 
occurred on Thursday, November 17th, at 6:00 p.m. The change in GPS week was 
expected at midnight on Saturday. The results looked fine and confirmed that the 
change had actually been made. The almanac data from week 776 was used, and the 
times from -55 to -49 hours corresponded to 113 to 1 19 hours of GPS week 775. It is 

noted that in the Novatel ephemeris file, the almanac entries read week 776, while the 

* 

navigation solution entries still read week 775. 
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Chapter 3 


Data Processing 

3.1 Introduction 

The goal is to observe and model the time and distance correlated behavior of 
DGPS errors. One such model is a first order Markov process, 


f 

MSE = Oqps 

V 


t 

T 


1-e 


y 


(3.1) 


where MSE is the mean square error, is the variance of the Markov process, x is 

the time correlation constant, and x is the distance correlation constant. A notional 

c 

plot of this model is shown in Figure 3.1. We will focus for the most part on Markov 
models because of their utility in navigation error analysis. The basic underlying phe- 
nomena which cause these errors may introduce correlation between the time and dis- 
tance behavior. An example would be a traveling ionospheric disturbance. Over an 


PREOECUNG 


35 



ensemble of cases, however, the correlation is expected to be zero. It will be so 
assumed for this study. 



Baseline Distance (km) 


Time Delay (seconds) 


Figure 3.1 Notional Mean Square Error Model for DGPS 


5000 


3.2 Inputs to the Processing 

The first step in the processing was to remove corrections for ionospheric delay. 
This eliminates the noise from the L 2 channel, reducing the noise level by a factor of 
^ . Differencing outputs from the two receivers will cancel the ionospheric errors to 
the extent that they are the same. This is a common practice in DGPS operation. Fur- 
ther characterization of this residual ionospheric error is one of the goals of the study. 
Removal of the ionospheric correction was actually done prior to the start of this study. 
[ 1 0] The format for the data is shown in Table 3.1. 


Table 3.1 GPS Data File Format 
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6 

7 

8 
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Time 

X 

Y 

Z 

# Sats 

Pml 

Pm2 

Pm3 

Pm4 










10 

11 

12 

13 

14 

15 

16 

17 

18 

Pm5 

Pm6 

GDOP 

PDOP 

HDOP 

VDOP 

TDOP 

Status 

Clock 
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The measurements from the two receivers A and B were edited by deleting mea- 
surements for which the time tags could not be matched. In addition, measurements 
were used only when the receivers indicated that the solution was good (status flag = 
1 ). 

Satellite switch times were computed in order to split the data into different sets 
that used the same four satellites for the entire set. This intentionally eliminates satel- 
lite switching as one on the sources of decorrelation. A negative aspect of this strategy 
is that we will be trying to estimate long term behavior from a number of short term 
data segments. This will be discussed in more detail later. The data segments were 
samples from what is presumed to be a zero mean process. The small sample size 
introduced a sample mean into the data. This mean was removed before determining 
the correlation between different data sets. 

Another strategy for getting as close to the underlying phenomenon as possible, 
i.e. errors in the pseudo-ranges, was to normalize the errors in the position solution 
components using the dilution of precision in each component, XDOP, YDOP, and 
ZDOR These values were computed by utilizing a GPS constellation simulation and 
inputting the measurement times for each set of data and the satellites tracked by each 
receiver. [6] The program makes use of the almanac file taken from the Novatel 
receivers. The DOP values were obtained from the diagonal of the (H^H) ^ matrix. 
The resulting DOPs were used to normalize the x, y, and z errors. 

For the rooftop, zero baseline data the GPS antenna was not placed at the surveyed 
marker. (The marker is in the “shadow” of a large vertical metal wall.) The placement 
was in an attempt to minimize multipath effects on the solution. Therefore, the truth 
value for the position (x, y, z) was defined to be the average position for all the data 
sets. This differed from the surveyed position by 5.6590 meters. 
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Ysurv meters, and + 18.7744 meters. These values were in the right 

direction and distance from the marker - as well as could be determined by eyesight. 
No effort was made to actually measure the distance, since it was noted that the survey 
marker position was itself only accurate to about a meter. [21] 

3.3 Definition of the Error and Difference in Errors 

This section defines the errors that were utilized in the autocorrelation processing. 

The quantity of interest is, finally, the difference in error between the time shifted 
errors from receiver A and/or B. 

^diff = (^crr (t + ” Xerr (0 ) /XDOP (3.2) 

The RSS difference is 

DIFFrss = /xJiff+yLf+4ff (3-3) 

The individual components of the error at each receiver are defined by; 

^err “ ^msd ^true (3-4) 

These individual errors are split into six sets where all the points in a set use the 
exact same four satellites. The RSS for each set was found using Equation (3.3). 

3.4 Auto and Cross Correlation 

The first strategy for characterizing the behavior of the differences was to examine 
the correlation in the rooftop data. The cross correlation and the autocorrelation of the 
rooftop, zero baseline, differences was computed. Since the two receivers were being 
driven from the same antenna, and the receivers were only contributing white noise, 
the cross and auto correlation were expected to both be measures of the same thing. 
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the behavior of the errors in pseudorange with time. This analysis was done for the 
component differences and the RSS difference for each of 6 data sets. 

The cross correlation between receiver A and B differences was found. The sample 
mean was subtracted from the data before the correlation algorithm was applied. Cor- 
relations were found for x, y, z, and the RSS differences for each set from 1 to 6 using 
Equation (3.2) and (3.3). 

The sigma and tau values were estimated for each set. The variance, a , is the cen- 
ter value, since the inputs were of zero mean. Tau is the value on the x axis that corre- 

2 

sponds to a y value of — . 

e 

A Markov process has a probability distribution that is dependent only on the 
value at one point immediately in the past. The process can be described by 
Equation (3.5), where w is white noise and p is the reciprocal time constant. [12] 


Jt 


+ P (t) X = w 


( 3 . 5 ) 


The autocorrelation of a zero mean Markov process is given by Equation (3.6). 
The correlation time, or 1/e point, is found at 1/p . [12] 




( 3 . 6 ) 


A single variance and time constant for the six set ensemble was found in two 
ways: a) by finding the correlation of the weighted average of the data in the six sets, 
and b) taking the weighted average of the variance and time constant found from the 
correlation of each of the six sets. 

The results of the cross-correlation process were disappointing for two reasons. 
The first is that a tremendous amount of data is required to get accurate results using 
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this technique. Brown and Hwang [7] present the limit on experimental autocorrela- 
tion as 


( 3 . 7 ) 

where i is the time constant, (x) is the experimentally determined autocorrelation 
function, T is the time length of the experimental record, and is the variance of the 
process. 

For a requirement of 10% accuracy on the estimation of the autocorrelation func- 
tion, the experimental autocorrelation should have a standard deviation less than 0.1 

2 « 4 VarrV (x)i 

times a , or 0.01 times a . Therefore, — 1 = 0.01 , and the limit becomes 

a 

^ ^ o.oip (^-8) 

for an accuracy of 10%. If the time constant is 1800 seconds, T = 360,000 sec- 
onds. In other words, 4.2 days of data is required. In general, the data required for 10% 
accuracy equals 200 times the time constant. A time constant of 250 seconds requires 
almost 14 hours of data for 10% accuracy. The development of the accuracy limit is 
presented in Appendix B. 

Figure 3.2 is a graphic expression of Equation (3.7). It shows the plot of amount of 
data vs. percent error in the autocorrelation function for a process with a time constant 
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of 1 800 seconds. Again, note that about four days worth of data is required in order to 


Amount of Data vs Percent Error for Time Constant of 1800 Seconds 



obtain about 10% error on the autocorrelation function. For our problem, we must 
have data where the satellites are all the same for each set of data. This limits the data 
sets to about one hour before a satellite switch occurs. The second reason for abandon- 
ing correlation as the analysis method was that it does not easily allow the resolution 
of multiple Markov processes. It will be found that there were indeed two separable 
Markov processes with, obviously, different time constants. 

3.5 Time Delay Error Analysis 

An alternative method for observing and characterizing the error inherent in using 
a time delayed DGPS correction was developed. It consisted of simply computing the 
time delay differences as a function of time delay and fitting those differences to a pre- 
sumed model. The time shift ranged from zero to a maximum such that 20 data points 
remained in common between the original set and the shifted set. Since the data inter- 
val was 10 seconds, this corresponds to an overlap of 200 seconds. Thus a 1500 sec- 
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ond data set will yield a 1300 second set of time shifted differences. The differences 
are the mean square, MS, of the components of the RMS differences at each time shift. 
After each shift, the shifted and unshifted errors were differenced. The differences 
(at each overlapping point) are functions of the total shift, x . (A shift by an index of 
one corresponds to a 10 second shift.) For the discrete case these differences are 

( 3 . 9 ) 

The mean and standard deviation at each shift, m, was found by summing over i 


X.5 


diff 


X * 


n. 


( 3 . 10 ) 




X diff = 


( 3 . 11 ) 


where the sum on i is over the number of overlapping points, at each shift, m. 
The mean square, MS, of the x, y, z standard deviations was computed 


= x^diff„ + y^diff„ + z^diff^ (3.12) 

where the mean square difference is expressed in terms of the time delay, x, corre- 
sponding to the shift, m. 
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Figure 3.3 shows this process in block diagram form. The result of the time delay 


error process is a curve of mean square error for differences vs. x . 



Figure 3.3 Creating Time Delay Errors 


3.6 Relation of Mean Square Differences to Markov Model 

Consider the outputs of receiver A and B: (x^, Ya* Za) and (xg, yg, Zg). Taking the 

RSS of each provides RSS^^ and RSSg as in Equation (3.10). The difference is DIFF = 
RSS^ - RSSg. The variance on the difference is given by 

^mFF ^ DIFF^-DTFF^ (3.13) 

If the difference is a zero mean process, then 

®mFF ' " 2RSS^RSSg + RSSg (3.14) 

Assume both A and B are zero mean processes, then RSS ^ , RSS^ = o„ , and 

A A 13 13 

RSS^RSSg = pcf^tJg . This leads to 

^DIFF ^ 
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^22 

Assuming A and B have similar variances, 0 ^ =: = 0 ^g resulting in 

Equation (3.16). 


Letting p = 


t X '^DIFF ” “ ^^AB ^ ^ 

X X 

Q 

e for a Markov process leads to the model function 


2 

^DIFF 





1 -e 


t 

X 


X 

X 


c 




which is equivalent to the mean square error if the bias is zero. 


( 3 . 16 ) 


( 3 . 17 ) 


3.7 Marquardt Method 

The Marquardt Method is a non-linear curve fitting routine that enables the fit of 
an exponential to a set of data. This is the means of obtaining the estimates for the 
DGPS model parameters. The inputs to the model include arrays of the independent 
and dependent data points, along with the standard deviations for each data point. The 
algorithm determines the standard deviation and correlation constant for the first order 
Markov fit, along with a covariance matrix that can be used to generate error bounds 
for the model. 

The process defines a chi-squared merit function and determines best-fit 
parameters by minimizing the merit function. For the nonlinear case, the minimization 
IS performed iteratively. The procedure is repeated until x effectively stops decreas- 
mg. The selected limit for the difference between consecutive values was 0. 1 . [22] 
The Marquardt method is described in more detail in Appendix A. The Marquardt 
method is used to determine independently the errors due to baseline distance and time 
delay. 


44 



3.8 Pseudorange Error Processing 

All previous processing was performed using the navigation solutions that are pro- 
vided by the GPS receiver. Due to the satellite switches, the data sets are limited in 
length. In an effort to obtain longer data sets, the possibility of obtaining pseudorange 
errors was investigated. The pseudorange from a single satellite could cross over sev- 
eral satellite switches and provide a source of continuous data for use in the process- 
ing. 

Equation (3.18) is used to obtain the pseudorange errors. 


p -p = H 
-meas -true 


r 

meas true 

At - At . , 

meas true 


where H is given by Equation (3.19). 


( 3 . 18 ) 


H = 


^1 ^ 
1 

1 

«3 ^ 
1 

«4 ^ 


( 3 . 19 ) 


The line of sight vector to each satellite is given by Ui^. The vector p contains 

-meas 

the measured pseudorange for each of the four satellites. The output GPS navigation 
solution vector is r . The receiver’s true position vector is r . The user clock 

meas * ‘true 

bias estimate is . The true user clock bias is At, . The line of sight vector to 

each satellite is obtained by using post-processed ephemeris data downloaded from the 
Internet. The ephemeris data is obtained at 15 minute intervals, and an interpolation 
routine called polint.m is used to find the true ephemeris for all data points [22]. 
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The problem with the pseudorange processing is that while we have a value for 
^tnie’ clock bias in Equation (3.18) is not known. If the error in the user clock 

bias, f^^ed to be zero, the pseudorange errors obtained will have 

jumps in them at each satellite switch. This is shown in the top half of Figure 3.4. If 
this jump is attributed to the change in At due to the satellite switch, a correction 

incBs 

term can be obtained which indeed allows the three continuing pseudoranges to be 
continuous. An example is shown in the bottom half of Figure 3.4. 

The actual clock error, however, is still unknown. We have simply found the 
change in its estimate due to the satellite switch. The difference between the measured 
and true clock error is a contributor to the pseudorange error, and cannot be ignored. It 
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may however be possible to estimate the contribution of this term statistically based on 
the TDOR This will be discussed in Section 5.2. 



3.9 Multipath 

As GPS signals are reflected off objects before entering the receiver antenna, a 
false pseudorange is created. This happens because the reflected signal has traveled a 
different path length than the direct signal. As noted in the previous section, cyclic 

behavior can be observed in the GPS position output which may be attributable to this 
“multipath” effect. 
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Oscillations are also observed in the time delay errors formed by shifting one 
receiver s data in time and differencing this with the other receiver’s output. Figure 3.5 
illustrates the point by showing the time delay error results from six data sets in the 
zero baseline case. Note the oscillations that occur in each set. 



After noting similar frequencies in each set, the average period was identified by 
taking the fast Fourier transform (FFl ) of each set shown. Figure 3.6 shows the FFT 
outputs from each of the six sets. The FFT output is very similar for all the sets with 
the tallest peaks occurring from 300 seconds to 600 seconds. Peaks off the scale were 
attributed to taking the FFT of a finite set of data. Peaks are expected to occur at 700 to 
1500 seconds because that’s how long the individual data segments lasted. The pre- 
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dominant periods which can reasonably be attributed to multipath are between 300 
seconds and 600 seconds. 



Period (sec) Period (sec) 

Figure 3.6 FFT of Time Delay Errors, Zero Baseline Data 


A scenario was developed which results in the multipath frequencies observed for 
the satellites in the zero baseline, 5-hour data set. The user selects a position offset 
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from the receiver for the location of the reflector. The true range and satellite position 
are used to determine the multipath frequency. The scenario is shown in Figure 3.7. 



Figure 3.7 Signal Multipath 


As shown, d] is the direct path from satellite to receiver. Two possible reflective 
paths, d 2 , were considered as examples. Note that d 2 includes the path before and after 
the reflection or the total from the satellite to the antenna. One was a reflection off a 
nearby building. The position offset used was 20 meters north and 20 meters up. The 
second multipath option is a reflection off the roof that enters under the receiver. This 
case is highly likely because the antenna was on a tripod and the gain under the 
receiver is fairly high. The offset used for this case was 1.9 meters down and 0.5 
meters north. 

The multipath frequency is given in Equation (3.20) and Equation (3.21). 


f(k) 


Ad, Ad2 
^At ^At 


( 3 . 20 ) 


f(k) 


d,(k+l) - d,(k) d2(k+l) -d2(k) 

0.19(t(k+l) -t(k)) ~ai9(t(k+l) -t(k)) 


( 3 . 21 ) 


where X is the wavelength of the signal (0.19 meters), t is the measurement time, and 
k is the index for the measurements. The results for satellite 26 for the first offset are 
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Multipath Frequency (Hz) 


shown in Figure 3.8. The period is between 90 and 100 seconds (1.5 to 1.7 minutes), 
and the frequency is 0.0095 to 0.011 Hz. 




GPS Time (hours) 

Figure 3.8 Sample Multipath Period for Offset of 20 Meters 
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Chapter 4 


Results 

4.1 Rooftop Experiment Results 

4.1.1 Navigation Solution for 5-Hour Data 

The first set of results is from processing the navigation solution for the 5-hour set 

of data with the zero baseline. Figure 4.1 shows the satellites that were used during the 
data collection. The y-axis is the satellite identification number. Note that at each time 
only four satellites are being used. The gap between 115 and 116.5 hours was noted 
earlier as the time when less than four satellites were tracked. 



The vertical lines mark each spot where there is a switch in at least one of the sat- 
ellites being tracked. The data was split into six different sets where the satellite set 
stayed constant over that interval. The small data segment towards the middle was not 
used. 


53 



Time delay differences were created using the processing described in Chapter 3. 
Figure 4.2 shows the resulting mean square vs. time delay plots for each of the six sets 
of data. The average is also plotted on this chart. 



Figure 4.3 shows the 1-sigma bounds on the time delay differences. These were 
computed by finding the standard deviation for the mean square error calculation as 
shown in Equation (4.1). [11] Note that the error bounds for the six sets overlap much 
of the time. 


^MS "* " 1) 


( 4 . 1 ) 
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2 

where is the mean square error (ordinate in graph), and n„, is the number of 
overlapping points. It can be seen that the error bounds increase as the number of over- 
lapping points decreases for increasing time delay. 



Figure 4.4 shows the weighted average from the six data sets, along with the stan- 
dard deviation on the average. These values are used as inputs to the Marquardt curve 
fitting routine. The weighted average and standard deviation were found using Equa- 
tion (4.2) and (4.3) where MSj is the mean square for each set and is the standard 
deviation for the mean square. 


mean 



( 4 . 2 ) 


55 



1 


c 


mean 


fl4 
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( 4 . 3 ) 



Figure 4.5 shows the Markov fit to the average from the 6 data sets. Equation (4.4) 
shows the model that was used in the Marquardt fitting routine. A sum of two Markov 
processes was used because this model provided a better fit to the data. The first pro- 
cess with the shorter time constant is attributed to multipath. The period falls in the 
sample range found in the multipath scenario as discussed in Chapter 3. All the zero 
baseline cases exhibit this short time constant behavior. The longer time constant pro- 
cess is intended to encompass all other errors. A model for this process was one of the 
goals of the study. 
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f 

Model = 

l-e^‘ 
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+ G 2 

1 - e 


^ J 
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( 4 . 4 ) 
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The error bounds on the Markov fit are shown with dotted lines. These are 
obtained by using the covariance provided by the Marquardt routine as shown in 
Equation (4.5). 


5e = . 


(4.5) 


Os 9s 

The one-sigma error on the model is 8e . The function f is the mean square error 
which is a function of the vector s . The vector s consists of both values for c and x . 
The partial derivatives of the function f with respect to the elements of s are output 
from the Marquardt routine. The routine also provides the covariance matrix E defined 
by Equation (4.6). 


T 

E = 5s5s (4.6) 

The noise level for the difference between receiver A and B is 0.05 meters^ and 
provides a bias for the navigation solution errors. The receiver noise has a mean of 
0.21 meters and a standard deviation of 0.23 meters. Squaring the standard deviation 
gives a bias for the mean square of 0.05 meters^. 

Figure 4.6 shows the residual for the Markov process fit. This is obtained by sub- 
tracting the fit from the data. The residual represents anything that has not been 
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Residual (meter^2) & Mean Square Error (melers'^2) 



uded in the model. Note the sinusoidal characteristic in the residual. This will be 
ressed later in the chapter. 



The model parameters for each of the 6 sets and the average are given in Table 4.1. 
The average values are obtained from a separate fit for the average curve, not just the 
average for the parameter values from the 6 data sets. The fits to the individual sets 
reveal the fact that the Markov model is not very sensitive to changes in the T 2 value. 
That is the reason for the large variations seen in this parameter. This insensitivity is 
related to the limit (variance) of the experimental correlation given in Chapter 3. 


Table 4.1 Results for Zero Baseline 


SET 



^2 

'^2 

(meters) 

(seconds) 

(meters) 

(seconds) 

1 

1.25 

272.5 

1.02 

285.9 

2 

0.68 

80.2 

5.74 

7643.3 

3 

0.85 

39.8 

26.4 

123,719 

4 

1.86 

204.9 

5.80 

11,412 

5 

1.47 

125.2 

3.92 

10,110 

6 

1.10 

57.6 

1.66 

506.6 

Average 

0.98 

103.9 

1.71 

499.0 


4,1.2 Navigation Solution for 24-Hour Data 

The second set of results is from processing the navigation solution for the 24-hour 

set of data with the zero baseline. Figure 4.7 shows the satellites that were used during 
the data collection. The y-axis is the satellite identification number. Note that at each 
time only four satellites are being used. The vertical lines mark each spot where there 
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ID Number 


is a switch in at least one of the satellites being tracked. Ten sections were selected for 
processing of the navigation solution. These are indicated by an ‘X’ above the data set. 



Time delay differences were created and Figure 4.8 shows the resulting mean 
square vs. time delay plots for each of the ten sets of data. The average is also plotted 
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on this chart. The average is the very thin line that runs through the middle of the 


curves. 



Noting that some of the curves in Figure 4.8 are lower than others, the time of day 
for each curve and the mean MSE were calculated. Figure 4.9 plots the mean of each 

of these sets of time delay differences versus time of day. The time of day is plotted 

* 

from 0 to 24 hours (1 days). Hour zero occurs at midnight. It might be said that there is 
a large mean error between 8:00 a.m. and 10:00 a.m. corresponding to a rise in total 
electron count in the morning, the smallest mean square differences occur between 
noon and 8:00 p.m. There is then another increase. This increase occurs significantly 
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after sunset. Any correlation between the mean square differences and sunrise and 
sunset is difficult to discern. 



Figure 4.10 shows the average from the ten data sets, along with the standard devi- 
ation on the average. These values are used as inputs to the Marquardt curve fitting 
routine. 



Figure 4.10 Average Navigation Solution Error, 24-Hour Data 
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Figure 4.11 shows the Markov fit to the average from the 10 data sets. The bias is 
zero because there is only one receiver for the 24-hour data. The result is a shift of 
receiver A vs. receiver A. The X 2 value of 3847.1 seconds equals 1.07 hours. 



Note that the error bounds on the raw data. Figure 4. 10, and on the Markov fit, Fig- 
ure 4.11, do not completely overlap. This is especially obvious after a few thousand 
seconds. This ‘inconsistency” reflects the fact that the Markov process is an incom- 
plete description of the pseudorange transmission delay differences. In fact, the errors 
and their differences are not stochastic, but are “unpredictable”. Even this is not 
exactly the case because, with enough effort, satellite ephemeris and transmission 
delay could be predicted at least in the short term. Characterizing these differences as 
Markov processes is merely convenient. As the analysis proceeds the character of the 
fit residuals will be observed. 


63 



Figure 4.12 shows the residual for the Markov fit on the 24-hour data. Again note 


the sinusoidal nature of the model error. 



Figure 4.13 shows the Markov fit for both the 5-hour and 24-hour data sets. The 
model parameters for these two sets are almost identical. The data was collected on 
different days and times of day. The results indicate a definite consistency between the 
data sets. Note that at the end of each data set the average error increases sharply. This 
is a result of the fact that the number of data points is decreasing at each shift. The 
smaller segments have only 20 points left at the end, or only 200 seconds of data. It 
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Mean Square Error (melers*2) 


seems reasonable that there would be trouble estimating a 1-hour time constant, with 
so few points remaining. 



The fit to each of the 10 navigation solution sets is given in Table 4.2. The large 
variance on the T 2 values is again apparent, and we will note once more that it is con- 
sistent with the variance in the correlation given in Chapter 3. 


Table 4.2 Results for Zero Baseline, 24-Hour Data 


SET 


'^1 

^2 

'^2 

(meters) 

(seconds) 

(meters) 

(seconds) 

1 

1,57 

89.8 

9.51 

38,358 

2 

1.81 

255.0 

12.54 

45,316 

3 

1.50 

208.6 

1.27 

214.0 

4 

1.07 

123.1 

41.4 

7.4x10^ 

5 

0.93 

105.1 

15.3 

46,620 

6 

1.98 

154.3 

1.36 

512.5 

7 

2.41 

296.7 

3.76 

22,579 

8 

0.69 

43.4 

13.1 

67,439 

9 

1.77 

252.0 

0.90 

261.6 


65 












































Table 4.2 Results for Zero Baseline, 24-Hour Data 


SET 

(meters) 

(seconds) 

^2 

(meters) 

(seconds) 

10 

1.54 

120.1 

32.1 

2.7x10^ 

Average 

1.47 

171.8 

2.77 

3847.1 


4.1.3 Pseudorange Errors, 5-Hour Data 

The third set of results is from processing the pseudorange errors for the 5-hour set 
of data with the zero baseline. Fourteen pseudorange sequences were formed from the 
data. The large gap in the data forced the splitting of pseudoranges on either side of the 
gap. Figure 4.14 shows the resulting mean square vs. time delay plots for each of the 
fourteen sets of data. The average is also plotted on this chart. 
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Mean Square Error (melers*2) _ Mean Square Error (meters^2j 


Figure 4.15 shows the average from the fourteen data sets, along with the standard 
deviation on the average. These values are used as inputs to the Marquardt curve fit- 
ting routine. 



Figure 4. 1 6 shows the Markov fit to the average from the 14 data sets. The first rise 


very close to what was seen in the navigation solution results. 
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Residual (meter^2) 


The residual for the Markov fit is shown in Figure 4.17. Once again note the sinu- 
soidal nature of the model error. 



Figure 4.17 Residual for Markov Fit, Pseudorange Errors 


4.1.4 Pseudorange Errors, 24-Hour Data 

The fourth set of results is from processing the pseudorange errors for the 24-hour 

set of data with the- zero baseline. Ten pseudorange sequences were formed from the 
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data. Figure 4.18 shows the resulting mean square vs. time delay plots for each of the 
ten sets of data. The average is also plotted on this chart. 
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Time Delay (seconds) 

Figure 4.18 Pseudorange Error, 24-Hour Data 


Figure 4.19 shows the average from the fourteen data sets, along with the standard 
deviation on the average. These values are used as inputs to the Marquardt curve fit- 
ting routine. 
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Figure 4.20 shows the Markov fit to the average from the 14 data sets. Again, The 
first rise is very close to what was seen in the navigation solution results. 



The residual for the Markov fit is shown in Figure 4.21. Once again note the sinu- 
soidal nature of the model error. 
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Mean Square Error (meters''2) 


The Markov fit for all cases is shown in Figure 4.22. All the curves are aligned in 
the first rise. 



4.1.5 Residual Errors 

All the residual errors shown have had sinusoidal characteristics. They are now 
examined more closely. Figure 4.23 plots the model residuals for all three of the previ- 
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ous cases. Not only are the sinusoidal characteristics evident in each, but the same fre- 
quency and phase is present. 



Figure 4.24 shows the FFT for each of the model residuals. The periods of the 
oscillations are very similar for all four cases. A primary period of about 400 seconds 
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is apparent. The next step is to see if the Markov model fits the data better if appropri- 
ate sinusoidal terms are added. 



100 


200 


300 


400 


500 


600 



A sinusoidal term using the period identified was added to the Markov model to 
attempt a better fit. The new fit is shown in Figure 4.25. This fit is much better than the 
Markov process alone. The model for the residual used is given in Equation (4.7). The 
short period of this oscillation once again points to multipath. Some particular, strong 
reflected signal common to all cases is a possible explanation 
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Mean Square Error (meters''2) 


Res = 0.35sin[27t[^j(t+ 150)1 (4.7) 



Figure 4.26 shows the residual for the new model. Note that the amplitude of the 
residual has been reduced to about a fourth of its former level. Note that sinusoidal 
characteristics can still be observed in the data. Further reduction using these periods 
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• Figure 4.28 shows the results from the West Coast data. 



The errors shown in Figure 4.27 and Figure 4.28 include the effects of multipath. 
An estimate of errors due to multipath was obtained using the same process that was 
used for the zero baseline data. That is, time shift differences were made, and a dual 
Markov model was fit to these differences. The shorter time constant was consistent 
with multipath effects. The variances were somewhat smaller than seen on the rooftop. 
Variances for each pair were obtained by combining the appropriate values. These 
variances are shown in Table 4.3. They were subtracted from the errors before pro- 
ceeding to find a combined time delay/distance model. 


Table 4.3 Multipath Variances for Each Baseline 


Separation (km) 

Locations 

®Liti (meters^) 

84 

Draper/Portsmouth 

0.22 

149 

Long Beach/San Diego 

0.92 

282 

Fallon/Stockton 

0.67 

496 

Fallon/Fort Irwin 

0.90 
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Mean Square Error (me1ers^2) 


Table 4.3 Multipath Variances for Each Baseline 


Separation (km) 

Locations 

oLiii (meters^) 

740 

Langley/Draper 

0.18 

756 

Fallon/San Diego 

0.49 

820 

Langley/Portsmouth 

0.10 


4.3 Combined Model 

The results have been presented for each of the three data sources: zero baseline. 
East Coast, and West Coast. Now the results from the navigation solution errors will 
be combined to create composite results. The 24-hour data for the zero baseline case 
will be used since it has a longer time base. 

4.3.1 Results From All Sets 

Figure 4.29 plots all baseline errors together. 
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Figure 4.30 plots the same information in 3-D format for a more visual interpreta- 
tion. 



Figure 4.30 3-D Error Plot, All Data Sets 


Figure 4.31 plots a cross section of the 3-D plot at zero time delay. This cross sec- 
tion will be used to fit the Markov process for the distance correlation. The last three 

points are from the East Coast baselines with a C/A-code receiver at one end (Willow 
♦ 

Run). The Draper data collection personnel indicated a lower confidence in the output 
from these receivers. In addition, the selective availability errors had to removed by 
Aerospace in post processing. By policy, these errors were not totally removed from 
the data, and some noise still remains. Given these facts, a strong case can be made for 
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Mean Square Error (meters^2) 


the removal of the C/A-code runs from the data sets. This would leave only P-code 
data leading to a higher confidence in the results. 
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Mean Square Error (melers^2) 


4.3.2 Results From High Confidence Sets 

To create a set of high confidence data, the three C/A-code sets at 838, 1022 and 

1052 km were removed. The resulting time delay errors are shown in Figure 4.32. 
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Figure 4.33 shows the 3-D plot for the high confidence data sets. Note that this sur- 


face is smoother than the one including the C/A-code results. 



Figure 4.33 3-D Error Plot, High Confidence Sets 


Figure 4.34 shows the 3-D plot for the high confidence data sets with the multipath 
removed from the zero baseline data. This plot will be compared with the final model. 



Figure 4.34 3-D Error Plot, No Multipath 
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Figure 4.35 plots the cross-sectional view of Figure 4.33 for the case with no time 
delay. This data is used as an input to the Marquardt method to fit the Markov model to 
the distance dimension. The 1 -sigma error bounds provided to the Marquardt method 
are also shown in the graph. The circles show the actual data points available for the 
Marquardt fit. Note that even with several baselines of collected data, there are still 
very few points in the distance dimension compared with the number of points in the 
time delay dimension. 



Figure 4.35 Zero Time Delay Results, High Confidence Sets 


Figure 4.36 shows the Markov fit for the distance dimension. The bias is zero 
because the 24-hour data was used at 0 km, and only one receiver was available for 
this data. Equation (4.8) shows the model that was used in the Marquardt fitting rou- 
tine for the distance dimension. 
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Figure 4.37 plots the residual for the distance dimension Markov fit. There is a 
slight sinusoidal characteristic to this plot, but not enough data points are available to 
merit trying to correct the model. 



Figure 4.38 shows the final resulting DGPS model obtained by combining the 
Markov fit to the time delay and baseline distance dimensions. Since the same under- 
laying phenomena cause the distance and time decorrelation, the variance of the pro- 
cess is common. The value of 3.73 meters^ is the weighted average of the values found 
for the time and distance models. Note that the raw data has a significant spread and 
this model will not encompass every case. Rather it is meant to provide an average 
error model that represents the typical behavior of the DGPS operation. The formula 
for the DGPS model is given by Equation (4.9). This plot can be compared with Figure 
4.34. Note that these results have been normalized by the DOP values. To determine 
the errors obtained in the navigation solution, the model must be multiplied by the 
appropriate DOP values for the given satellite geometry. 
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MS = 3.73 


t 


1 - e 


3847.1 122.8 


( 4 . 9 ) 



Figure 4.38 DGPS Error Model 


The data used to determine the time constant, t , comes exclusively from the “roof- 
top experiments. This is because the data from the East and West Coast experiments 
was not of sufficient duration to contribute much to the determination of T . 
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Chapter 5 


Conclusions and Suggestions for Future Work 


This thesis analyzed GPS data collected by Draper Laboratory to determine DGPS 
error characteristics. It was found that a first order Markov process could be used as an 
appropriate model for the errors in differential GPS. Since the errors are correlated in 
both time and distance, the Markov model has a process variance, a time correlation 
constant, and a distance correlation constant. Analysis of both navigation solution 
errors and pseudorange bias errors provided a wide range of data sets for processing. 
The Marquardt non-linear curve fitting routine was used to obtain the parameters for 
the Markov model for each set of data. 


5.1 Conciusions 

The resulting DGPS model is given in Equation (5.1). 




MS = 3.73 


1 -e 


t X 

3847.1 ~ 122.8 


meters^ 


V 


(5.1) 


where MS is the mean square error in the pseudorange biases. 


This model is fairly representative of average navigation solution errors. It should 
be noted that there are significant deviations from nominal performance. (Refer to Fig- 
ure 4.8 for example.) The analyst might wish to choose a larger value for the variance 
to describe “worst case” performance. The rationale for using a Markov fit is that 
many physical processes tend to exhibit Markov-like behavior and that it is a conve- 
nient model for linearized analysis. 


PiECQVHG 
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The error model is consistent with the differential GPS error budget given in Chap- 
ter 1. The DGPS budget had a MS error of 1.59 meters^, (1.26^). The Markov model 
detenmned herein has a MS error of 3.73 meters^. The variance at 92.6 km separation 
and zero time delay is about 1,98 meters^. Fifty nautical miles (92.6 kilometers) is the 
distance assumed in the referenced error budget. Multipath effects could easily 
account for the difference. The referenced budget is only the specification for the sys- 
tem and not necessarily indicative of actual performance. 

The raw data has a significant spread to it depending on the time of day, geometry, 
multipath, weather conditions, and possibly other factors. This model does not encom- 
pass every facet of the error sources. In particular we have seen the importance of mul- 
tipath error modeling. An attempt was also made to correlate the magnitude of the 
time delay errors to the time of day. Few data points were available, but the results 
indicate the lowest values occur during the afternoon and early evening. This would 
mean that the output from the two receivers is more highly correlated at these times. 

The model was developed by analyzing the errors in the navigation solution output 
of the GPS receivers. The individual pseudorange errors that contribute to the position 
error were also examined. A technique was developed to attribute changes in all pseu- 
dorange errors at satellite switch to changes in the user clock bias error estimate. This 
indeed caused the pseudorange error for each satellite to be continuous and not experi- 
ence jumps when a satellite switch occurred. The results from this analysis were con- 
sistent with those obtained from the navigation solution analysis for the zero baseline 
data. While this is encouraging, and lends more credibility to the determination of the 
time constant, the lack of a measure for true GPS time means the value for pseudor- 
ange error determined in this manner will be consistently low. 
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The cross-correlation technique was ineffective at determining the Markov fit 
parameters. This is due to the requirement for large continuous data sets for cross-cor- 
relation accuracy. DGPS data sets are limited to about one hour due to satellite 
switches. Alternatively, the Marquardt curve fitting routine was used to fit the time 
delay error to a Markov process. This routine also allowed fitting of two independent 
processes. 

This thesis effort has revealed many interesting aspects of the GPS receiver opera- 
tion including a close look at P-code data from two receivers attached to the same 
antenna. Of biggest importance may be the dramatic change in amplitude and fre- 
quency of oscillations in the navigation solution at each point where a satellite switch 
occurs. Given that multipath is the cause of such behavior, shielding or other counter- 
measures should be investigated on future experiments or in operational scenarios. 

5.2 Future Work 

As just mentioned, multipath should be treated more carefully in future experi- 
ments. It should be eliminated in order to examine other error sources and should be 
introduced in a controlled fashion if desired. 

Internet GPS data may be able to provide longer sets of data where the restriction 
of common satellites in each set is no longer a requirement. In order to obtain this, the 
best sites with accurately surveyed positions and hydrogen maser clocks would pro- 
vide the best opportunity for such data. This data would provide a check against the 
work done in this thesis and provide more insight into the correlation constants 
obtained. The best Internet address to use for these sources is toba.ucsd.edu. It is not 
clear, however, how the errors such as selective availability can be completely sepa- 
rated from the other physical errors of interest. The Internet receivers are not 
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encrypted, so SA will be present. Since the position is known very accurately, the 
errors in the navigation solution will be known very accurately, but the contribution 
from each error source can only be estimated. 

Even if true GPS time is unavailable at the receiver, a comparison of Time Dilution 
of Precision (TDOP) values with the PDOP values could be utilized to obtain a statis- 
tical estimate of the receiver clock errors. This would allow inclusion of these errors 
into the estimate of the pseudorange errors combined with the adjustment of clock bias 
devised herein. This would allow longer data sets to be analyzed. Short segments of 
data were one of the primary limitations for this thesis effort. 
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Appendix A 


Marquardt Method 

The Marquardt method or Levenberg-Marquardt method provides a way to per- 
form non-linear curve fitting for data. The process defines a chi-squared merit 
function and determines best-fit parameters by minimizing the merit function. For the 

nonlinear case, the minimization is performed iteratively. The procedure is repeated 
2 

until X effectively stops decreasing. Numerical Recipes [22] explains the Marquardt 
method in detail. The main points will be discussed here. 

The model to be fitted is 


y = y(Xj;a) 


(A.2) 


where a is the vector of parameters in the model that the Marquardt method must find 


In this 


is case, a = [a and y = a 


1-e ^ 


2,. 


. The X merit function is 


(5) = x 


i= 1 


Xj-y (Xj;a) 


a \-\2 


(A.3) 


The goal of the Marquardt method is to find a„,in that minimizes x^ • An initial guess 
called a cur is required to get things started. If the guess for a is fairly good, the second 
derivative matrix (Hessian matrix) of the x function is used to jump immediately to 
the a for that guess using 


n = + • [-Vx^(a^j] 


(A.4) 


where D is the Hessian matrix. 
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However, if the guess is poor, the inverse-Hessian method cannot be used. In this 
case, the steepest descent method is used to take a step down the gradient 


^next = ^cur - constant xVx“(a^J (A.5) 

where the constant is small enough not to exhaust the downhill direction. 

2 

The gradient of % with respect to the parameters a , which will be zero at the 
minimum, has components 


3a,^ ^ 

i= 1 


j n 

2 (Xi;a) k = 1, 2, ..., M 


(A.6) 


Taking an additional partial derivative gives 


N 


1 s 1 1 K I ^ 

The following terms are defined 


(A.7) 


3 =-i^ 

23a. 


ttu, S - 


2 0 


23ajj3aj 


(A.8) 


Now [a] — -D in Equation (A.4) which can be rewritten 


M 

^ttkiSaj = 3k (A.9) 

i = I 

This set of linear equations can be solved for the increments 5aj that are added to the 
current approximation to get the next approximation. 

The steepest descent formula translates to 

6aj = constant x3| (A.iO) 
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The Marquardt method allows the smooth variation between the steepest descent 
method and the inverse-Hessian method as the approximation goes from poor to good. 
The two methods are combined into one equation 

M 

2^aki5a, = (A.11) 

i = 1 

where a is given by 

a'j s ( 1 + ?L) 

, (A.12) 

(J^k) 

and X is referred to by the author as a non-dimensional fudge factor. If X is large, the 

approximation is poor, and the steepest descent method is used. If X is near zero, the 

approximation is good, and the inverse-Hessian method is used. 

The Marquardt recipe consists of the following: 

2 

• Compute % (a) 

• Pick a modest value for X, say ^ = 0.001. 

• (1) Solve Equation (A. 11) for 5a and ev 2 iluate (a + 5a) 

• If X (a + 5a) > x^ (a) , increase A. by a factor of 10 and return to (1) 

• If X (a + 5a) < X (a) , decrease A- by a factor of 10, update the trial solution 
a <— a + 5a , and return to (1). 

Once the acceptable minimum is found (change in % of less than 0.1), A. is set to 
zero and the covariance is found 


[C] = [a]-' (A. 13 ) 

The partial derivatives of the model with respect to the individual parameters in a 
are used by the Marquardt method. For the Markov exponential fit, the partials are 


^ = 
Oa 


/ 

2c 



V ; 


(A.14) 
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(A. 1 5) 



Appendix B 


Autocorrelation Limits 

The limits on the experimental autocorrelation function were found in [7]. The ini- 
tial assumption proven in [3] is that the variance of an experimentally determined 
autocorrelation function satisfies the inequality 

oe 

VarV,(x) <|jR^(T)dT (B.i) 

0 

where T = time length of the experimental record 

(t) = autocorrelation function of the Gaussian process under consideration 
VJ'C) = time average of X-r(t)XT (t + T) where X^-Ct) is the finite-length sam- 
ple of X (t) (that is, (t) is the experimentally determined autocorrelation function 
based on a finite record length) 

The experimental autocorrelation is formed by 

= f^/r (‘ + -t) dt (B.2) 

It is noted here that the Matlab function xcorr does not divide by — 5 — , so the user 
must be careful when working with that function. Dividing properly will result in 
Equation (B.2). 

(t) can be shown to be an unbiased estimator of (r) [7]. Assuming that the 
process X(t) is a Gauss-Markov process, then the autocorrelation function is 

Rx (T) = c^e (B.3) 

Note that P = . This autocorrelation can be substituted in Equation (B.I) to give 
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(B.4) 


which is the limit that is used for the accuracy on an experimental autocorrelation 
function in Chapter 3. 
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